#include "functionSpace.h"

using namespace std;

static double Lagrange_Tri_refNodesPos_1    [3][2]  = {{0,0},{1,0},{0,1}};   
static int    Lagrange_Tri_faceNodesIndex_1 [3][2]  = {{0,1},{1,2},{2,0}};
static double Lagrange_Tri_refNodesPos_2    [6][2]  = {{0,0},{1,0},{0,1},{0.5,0},{0.5,0.5},{0,0.5}};  
static int    Lagrange_Tri_faceNodesIndex_2 [3][3]  = {{0,3,1},{1,4,2},{2,5,0}};
static double Lagrange_Tri_refNodesPos_3    [10][2] = {{0,0},{1,0},{0,1},{1./3.,0},{2./3.,0},{2./3.,1./3.},{1./3.,2./3.},{0,2./3.},{0,1./3.},{1./3.,1./3.}};   
static int    Lagrange_Tri_faceNodesIndex_3 [3][4]  = {{0,3,4,1},{1,5,6,2},{2,7,8,0}};

static double Tri_QP_1 [1][2]  = {{0.333333333333,0.333333333333}};
static double Tri_w_1  [1]     = {1};
static double Tri_QP_2 [3][2]  = {{0.166666666667,0.166666666667},{0.166666666667,0.666666666667},{0.666666666667,0.166666666667}};
static double Tri_w_2  [3]     = {0.333333333333,0.333333333333,0.333333333333};
static double Tri_QP_3 [4][2]  = {{0.333333333333,0.333333333333},{0.2,0.2},{0.2,0.6},{0.6,0.2}};
static double Tri_w_3  [4]     = {-0.5625,0.520833333333,0.520833333333,0.520833333333};
static double Tri_QP_4 [6][2]  = {{0.445948490916,0.445948490916},{0.445948490916,0.108103018168},{0.108103018168,0.445948490916},{0.0915762135098,0.0915762135098},{0.0915762135098 ,0.81684757298},{0.81684757298,0.0915762135098}};
static double Tri_w_4  [6]     = {0.223381589678,0.223381589678,0.223381589678,0.109951743655,0.109951743655,0.109951743655};
static double Tri_QP_5 [7][2]  = {{0.333333333333,0.333333333333},{0.470142064105,0.470142064105},{0.470142064105,0.0597158717898},{0.0597158717898,0.470142064105},{0.101286507323,0.101286507323},{0.101286507323,0.797426985353},{0.797426985353,0.101286507323}};
static double Tri_w_5  [7]     = {0.225,0.132394152789,0.132394152789,0.132394152789,0.125939180545,0.125939180545,0.125939180545};
static double Tri_QP_6 [12][2] = {{0.249286745171,0.249286745171},{0.249286745171,0.501426509658},{0.501426509658,0.249286745171},{0.0630890144915,0.0630890144915},{0.0630890144915,0.873821971017},{0.873821971017,0.0630890144915},{0.310352451034,0.636502499121},{0.636502499121,0.0531450498448},{0.0531450498448,0.310352451034},{0.310352451034,0.0531450498448},{0.636502499121,0.310352451034},{0.0531450498448,0.636502499121}};
static double Tri_w_6  [12]    = {0.116786275726,0.116786275726,0.116786275726,0.0508449063702,0.0508449063702,0.0508449063702,0.0828510756184,0.0828510756184,0.0828510756184,0.0828510756184,0.0828510756184,0.0828510756184};
static double Tri_QP_7 [13][2]  = {{0.333333333333,0.333333333333},{0.260345966079,0.260345966079},{0.260345966079,0.479308067842},{0.479308067842,0.260345966079},{0.0651301029022 ,0.0651301029022},{0.0651301029022,0.869739794196},{0.869739794196,0.0651301029022},{0.312865496005,0.63844418857},{0.63844418857,0.0486903154253},{0.0486903154253,0.312865496005},{0.312865496005,0.0486903154253},{0.63844418857,0.312865496005},{0.0486903154253,0.63844418857}};
static double Tri_w_7  [13]     = {-0.149570044468,0.175615257433,0.175615257433,0.175615257433,0.0533472356088,0.0533472356088,0.0533472356088,0.0771137608903,0.0771137608903,0.0771137608903,0.0771137608903,0.0771137608903,0.0771137608903};
static double Tri_QP_8 [16][2]  = {{0.333333333333,0.333333333333},{0.459292588293,0.459292588293},{0.459292588293,0.0814148234146},{0.0814148234146,0.459292588293},{0.170569307752 ,0.170569307752},{0.170569307752,0.658861384496},{0.658861384496,0.170569307752},{0.050547228317,0.050547228317},{0.050547228317,0.898905543366},{0.898905543366,0.050547228317},{0.263112829635,0.728492392955},{0.728492392955,0.00839477740996},{0.00839477740996,0.263112829635},{0.263112829635,0.00839477740996},{0.728492392955,0.263112829635},{0.00839477740996,0.728492392955}};
static double Tri_w_8  [16]     = {0.144315607678,0.0950916342673,0.0950916342673,0.0950916342673,0.103217370535,0.103217370535,0.103217370535,0.0324584976232,0.0324584976232,0.0324584976232,0.0272303141744,0.0272303141744,0.0272303141744,0.0272303141744,0.0272303141744,0.0272303141744};
static double Tri_QP_9 [19][2]  = {{0.333333333333,0.333333333333},{0.489682519199,0.489682519199},{0.489682519199,0.0206349616025},{0.0206349616025,0.489682519199},{0.437089591493 ,0.437089591493},{0.437089591493,0.125820817014},{0.125820817014,0.437089591493},{0.188203535619,0.188203535619},{0.188203535619,0.623592928762},{0.623592928762,0.188203535619},{0.0447295133945,0.0447295133945},{0.0447295133945,0.910540973211},{0.910540973211,0.0447295133945},{0.221962989161,0.741198598784},{0.741198598784,0.0368384120547},{0.0368384120547,0.221962989161},{0.221962989161,0.0368384120547},{0.741198598784,0.221962989161},{0.0368384120547,0.741198598784}};
static double Tri_w_9  [19]     = {0.0971357962828,0.0313347002271,0.0313347002271,0.0313347002271,0.0778275410048,0.0778275410048,0.0778275410048,0.0796477389272,0.0796477389272,0.0796477389272,0.0255776756587,0.0255776756587,0.0255776756587,0.0432835393773,0.0432835393773,0.0432835393773,0.0432835393773,0.0432835393773,0.0432835393773};
static double Tri_QP_10 [25][2]  = {{0.333333333333,0.333333333333},{0.485577633384,0.485577633384},{0.485577633384,0.0288447332327},{0.0288447332327,0.485577633384},{0.109481575485,0.109481575485},{0.109481575485,0.78103684903},{0.78103684903,0.109481575485},{0.307939838764,0.550352941821},{0.550352941821,0.141707219415},{0.141707219415,0.307939838764},{0.307939838764,0.141707219415},{0.550352941821,0.307939838764},{0.141707219415,0.550352941821},{0.24667256064,0.728323904597},{0.728323904597,0.0250035347627},{0.0250035347627,0.24667256064},{0.24667256064,0.0250035347627},{0.728323904597,0.24667256064},{0.0250035347627,0.728323904597},{0.0668032510122,0.923655933587},{0.923655933587,0.0095408154003},{0.0095408154003,0.0668032510122},{0.0668032510122,0.0095408154003},{0.923655933587,0.0668032510122},{0.0095408154003,0.923655933587}};
static double Tri_w_10  [25]     = {0.0908179903828,0.0367259577565,0.0367259577565,0.0367259577565,0.0453210594355,0.0453210594355,0.0453210594355,0.0727579168454,0.0727579168454,0.0727579168454,0.0727579168454,0.0727579168454,0.0727579168454,0.0283272425311,0.0283272425311,0.0283272425311,0.0283272425311,0.0283272425311,0.0283272425311,0.00942166696373,0.00942166696373,0.00942166696373,0.00942166696373,0.00942166696373,0.00942166696373};
static double Tri_QP_11 [27][2]  = {{0.534611048271,0.534611048271},{0.534611048271,-0.0692220965415},{-0.0692220965415,0.534611048271},{0.398969302966,0.398969302966},{0.398969302966,0.202061394068},{0.202061394068,0.398969302966},{0.203309900431,0.203309900431},{0.203309900431,0.593380199137},{0.593380199137,0.203309900431},{0.119350912283,0.119350912283},{0.119350912283,0.761298175435},{0.761298175435,0.119350912283},{0.0323649481113,0.0323649481113},{0.0323649481113,0.935270103777},{0.935270103777,0.0323649481113},{0.356620648261,0.593201213428},{0.593201213428,0.0501781383105},{0.0501781383105,0.356620648261},{0.356620648261,0.0501781383105},{0.593201213428,0.356620648261},{0.0501781383105,0.593201213428},{0.171488980304,0.80748900316},{0.80748900316,0.0210220165362},{0.0210220165362,0.171488980304},{0.171488980304,0.0210220165362},{0.80748900316,0.171488980304},{0.0210220165362,0.80748900316}};
static double Tri_w_11  [27]     = {0.000927006328961,0.000927006328961,0.000927006328961,0.0771495349148,0.0771495349148,0.0771495349148,0.0593229773808,0.0593229773808,0.0593229773808,0.0361845405034,0.0361845405034,0.0361845405034,0.0136597310027,0.0136597310027,0.0136597310027,0.0523371119622,0.0523371119622,0.0523371119622,0.0523371119622,0.0523371119622,0.0523371119622,0.0207076596391,0.0207076596391,0.0207076596391,0.0207076596391,0.0207076596391,0.0207076596391};

static double Tri_coef_1 [3][3]   = {{1,-1,-1},{0,1,0},{0,0,1}};
static double Tri_coef_2 [6][6]   = {{1,-3,-3,2,2,4},{0,-1,0,2,0,0},{0,0,-1,0,2,0},{0,4,0,-4,0,-4},{0,0,0,0,0,4},{0,0,4,0,-4,-4}};
static double Tri_coef_3 [10][10] = {{1,-5.5,-5.5,9,9,18,-4.5,-4.5,-13.5,-13.5},{0,1,0,-4.5,0,0,4.5,0,0,0},{0,0,1,0,-4.5,0,0,4.5,0,0},{0,9,0,-22.5,0,-22.5,13.5,0,27,13.5},{0,-4.5,0,18,0,4.5,-13.5,0,-13.5,0},{0,0,0,0,0,-4.5,0,0,13.5,0},{0,0,0,0,0,-4.5,0,0,0,13.5},{0,0,-4.5,0,18,4.5,0,-13.5,0,-13.5},{0,0,9,0,-22.5,-22.5,0,13.5,13.5,27},{0,0,0,0,0,27,0,0,-27,-27}};

static double Edge_coef_1 [2][2]  = {{0.5,-0.5},{0.5,0.5}};
static double Edge_coef_2 [3][3]  = {{0,-0.5,0.5},{1,0,-1},{0,0.5,0.5}};
static double Edge_coef_3 [4][4]  = {{-0.0625,0.0625,0.5625,-0.5625},{0.5625,-1.6875,-0.5625,1.6875},{0.5625,1.6875,-0.5625,-1.6875},{-0.0625,-0.0625,0.5625,0.5625}};

static double Edge_QP_1[1] = {0.000000000000000e+00};
static double Edge_QP_2[2] = {-5.773502691896257e-01,5.773502691896257e-01};
static double Edge_QP_3[3] = {-7.745966692414834e-01,0.000000000000000e+00,7.745966692414834e-01};
static double Edge_QP_4[4] = {-8.611363115940526e-01,-3.399810435848563e-01,3.399810435848563e-01,8.611363115940526e-01};
static double Edge_QP_5[5] = {-9.061798459386640e-01,-5.384693101056831e-01,0.000000000000000e+00,5.384693101056831e-01,9.061798459386640e-01};
static double Edge_QP_6[6] = {-9.324695142031521e-01,-6.612093864662646e-01,-2.386191860831969e-01,2.386191860831969e-01,6.612093864662646e-01,9.324695142031521e-01};
static double Edge_QP_7[7] = {-9.491079123427585e-01,-7.415311855993945e-01,-4.058451513773972e-01,0.000000000000000e+00,4.058451513773972e-01,7.415311855993945e-01,9.491079123427585e-01};

static double Edge_w_1[1] = {2.000000000000000e+00};
static double Edge_w_2[2] = {1.000000000000000e+00,1.000000000000000e+00};
static double Edge_w_3[3] = {5.555555555555552e-01,8.888888888888888e-01,5.555555555555552e-01};
static double Edge_w_4[4] = {3.478548451374537e-01,6.521451548625464e-01,6.521451548625464e-01,3.478548451374537e-01};
static double Edge_w_5[5] = {2.369268850561890e-01,4.786286704993665e-01,5.688888888888889e-01,4.786286704993665e-01,2.369268850561890e-01};
static double Edge_w_6[6] = {1.713244923791705e-01,3.607615730481386e-01,4.679139345726913e-01,4.679139345726913e-01,3.607615730481386e-01,1.713244923791705e-01};
static double Edge_w_7[7] = {1.294849661688697e-01,2.797053914892767e-01,3.818300505051190e-01,4.179591836734694e-01,3.818300505051190e-01,2.797053914892767e-01,1.294849661688697e-01};

static double pp[10][2] = {{0,0},{1,0},{0,1},{2,0},{0,2},{1,1},{3,0},{0,3},{2,1},{1,2}};


//------------------------------------------------------------------------------
LagrangeFctSpaceForTriangle::LagrangeFctSpaceForTriangle(int order){

  this->order  = order;
  this->Qorder = 2*order+1;
  
  setRefNodes();
  setQuadratureRules();
  setMatrices();
    
}
//------------------------------------------------------------------------------
void LagrangeFctSpaceForTriangle::setRefNodes(){
  
  switch (this->order) {
    case 1:{
      // 2         | \ 
      // |\        f3 f2
      // 0-1       |-f1-
      nbENP          = 3;
      nbFNP          = 2;    
      ERefNodesPos = DGMatrix(nbENP,2);
      for(int i=0;i<nbENP;i++) for(int j=0;j<2;j++)      ERefNodesPos(i,j) = Lagrange_Tri_refNodesPos_1 [i][j];
      DGIntMatrix* tmp = new DGIntMatrix(3,nbFNP); FRefNodesIndex = *tmp;
      for(int i=0;i<3;i++) for(int j=0;j<nbFNP;j++)    FRefNodesIndex(i,j) = Lagrange_Tri_faceNodesIndex_1 [i][j];
      ERefNodesCoef = DGMatrix(nbENP,nbENP);
      for(int i=0;i<nbENP;i++) for(int j=0;j<nbENP;j++) ERefNodesCoef(i,j) = Tri_coef_1 [i][j]; 
      FRefNodesCoef = DGMatrix(nbFNP,nbFNP);
      for(int i=0;i<nbFNP;i++) for(int j=0;j<nbFNP;j++) FRefNodesCoef(i,j) = Edge_coef_1 [i][j];
      break;}
    case 2:{
      // 2
      // |\
      // 5 4
      // |  \
      // 0-3-1
      nbENP          = 6;
      nbFNP          = 3;      
      ERefNodesPos = DGMatrix(nbENP,2);
      for(int i=0;i<nbENP;i++) for(int j=0;j<2;j++)      ERefNodesPos(i,j) = Lagrange_Tri_refNodesPos_2 [i][j];
      DGIntMatrix* tmp = new DGIntMatrix(3,nbFNP); FRefNodesIndex = *tmp;
      for(int i=0;i<3;i++) for(int j=0;j<nbFNP;j++)    FRefNodesIndex(i,j) = Lagrange_Tri_faceNodesIndex_2 [i][j];
      ERefNodesCoef = DGMatrix(nbENP,nbENP);
      for(int i=0;i<nbENP;i++) for(int j=0;j<nbENP;j++) ERefNodesCoef(i,j) = Tri_coef_2 [i][j]; 
      FRefNodesCoef = DGMatrix(nbFNP,nbFNP);
      for(int i=0;i<nbFNP;i++) for(int j=0;j<nbFNP;j++) FRefNodesCoef(i,j) = Edge_coef_2 [i][j];
      break;}
    case 3:{
      // 2
      // |\
      // 7 6
      // |  \
      // 8 9 5
      // |    \
      // 0-3-4-1
      nbENP          = 10;
      nbFNP          = 4;      
      ERefNodesPos = DGMatrix(nbENP,2);
      for(int i=0;i<nbENP;i++) for(int j=0;j<2;j++)      ERefNodesPos(i,j) = Lagrange_Tri_refNodesPos_3 [i][j];
      DGIntMatrix* tmp = new DGIntMatrix(3,nbFNP); FRefNodesIndex = *tmp;
      for(int i=0;i<3;i++) for(int j=0;j<nbFNP;j++)    FRefNodesIndex(i,j) = Lagrange_Tri_faceNodesIndex_3 [i][j];
      ERefNodesCoef = DGMatrix(nbENP,nbENP);
      for(int i=0;i<nbENP;i++) for(int j=0;j<nbENP;j++) ERefNodesCoef(i,j) = Tri_coef_3 [i][j]; 
      FRefNodesCoef = DGMatrix(nbFNP,nbFNP);
      for(int i=0;i<nbFNP;i++) for(int j=0;j<nbFNP;j++) FRefNodesCoef(i,j) = Edge_coef_3 [i][j];
      break;}
    default:
      break;
  }
  
  Epp = DGMatrix(nbENP,2);
  for(int i=0;i<nbENP;i++) for(int j=0;j<2;j++) Epp(i,j) = pp[i][j];
  
}
//------------------------------------------------------------------------------
void LagrangeFctSpaceForTriangle::setQuadratureRules(){

  switch (this->Qorder) {
    case 1:
      nbEQP  = 1;
      nbFQP  = 1;
      EQP = DGMatrix(nbEQP,2);
      EW  = DGVector(nbEQP);
      FQP = DGVector(nbFQP);
      FW  = DGVector(nbFQP);
      for(int i=0;i<nbEQP;i++){
                             EW(i)    = Tri_w_1   [i]; 
        for(int j=0;j<2;j++) EQP(i,j) = Tri_QP_1  [i][j];
        if(i<nbFQP)         {FQP(i)   = Edge_QP_1 [i];
                             FW(i)    = Edge_w_1  [i];}
      }
      break;
    case 2:
      nbEQP  = 3;
      nbFQP  = 2;
      EQP = DGMatrix(nbEQP,2);
      EW  = DGVector(nbEQP);
      FQP = DGVector(nbFQP);
      FW  = DGVector(nbFQP);
      for(int i=0;i<nbEQP;i++){
                             EW(i)    = Tri_w_2   [i]; 
        for(int j=0;j<2;j++) EQP(i,j) = Tri_QP_2  [i][j];
        if(i<nbFQP)         {FQP(i)   = Edge_QP_2 [i];
                             FW(i)    = Edge_w_2  [i];}
      }
      break;
    case 3:
      nbEQP  = 4;
      nbFQP  = 3;
      EQP = DGMatrix(nbEQP,2);
      EW  = DGVector(nbEQP);
      FQP = DGVector(nbFQP);
      FW  = DGVector(nbFQP);
      for(int i=0;i<nbEQP;i++){
                             EW(i)    = Tri_w_3   [i]; 
        for(int j=0;j<2;j++) EQP(i,j) = Tri_QP_3  [i][j];
        if(i<nbFQP)         {FQP(i)   = Edge_QP_3 [i];
                             FW(i)    = Edge_w_3  [i];}
      }
      break;
    case 4:
      nbEQP  = 6;
      nbFQP  = 4;
      EQP = DGMatrix(nbEQP,2);
      EW  = DGVector(nbEQP);
      FQP = DGVector(nbFQP);
      FW  = DGVector(nbFQP);
      for(int i=0;i<nbEQP;i++){
                             EW(i)    = Tri_w_4   [i]; 
        for(int j=0;j<2;j++) EQP(i,j) = Tri_QP_4  [i][j];
        if(i<nbFQP)         {FQP(i)   = Edge_QP_4 [i];
                             FW(i)    = Edge_w_4  [i];}
      }
      break;
    case 5:
      nbEQP  = 7;
      nbFQP  = 5;
      EQP = DGMatrix(nbEQP,2);
      EW  = DGVector(nbEQP);
      FQP = DGVector(nbFQP);
      FW  = DGVector(nbFQP);
      for(int i=0;i<nbEQP;i++){
                             EW(i)    = Tri_w_5   [i]; 
        for(int j=0;j<2;j++) EQP(i,j) = Tri_QP_5  [i][j];
        if(i<nbFQP)         {FQP(i)   = Edge_QP_5 [i];
                             FW(i)    = Edge_w_5  [i];}
      }
      break;
    case 6:
      nbEQP  = 12;
      nbFQP  = 6;
      EQP = DGMatrix(nbEQP,2);
      EW  = DGVector(nbEQP);
      FQP = DGVector(nbFQP);
      FW  = DGVector(nbFQP);
      for(int i=0;i<nbEQP;i++){
                             EW(i)    = Tri_w_6   [i]; 
        for(int j=0;j<2;j++) EQP(i,j) = Tri_QP_6  [i][j];
        if(i<nbFQP)         {FQP(i)   = Edge_QP_6 [i];
                             FW(i)    = Edge_w_6  [i];}
      }
      break;
    case 7:
      nbEQP  = 13;
      nbFQP  = 7;
      EQP = DGMatrix(nbEQP,2);
      EW  = DGVector(nbEQP);
      FQP = DGVector(nbFQP);
      FW  = DGVector(nbFQP);
      for(int i=0;i<nbEQP;i++){
                             EW(i)    = Tri_w_7   [i]; 
        for(int j=0;j<2;j++) EQP(i,j) = Tri_QP_7  [i][j];
        if(i<nbFQP)         {FQP(i)   = Edge_QP_7 [i];
                             FW(i)    = Edge_w_7  [i];}
      }
      break;
    default:
      break;
  }
  
}
//------------------------------------------------------------------------------
double LagrangeFctSpaceForTriangle::getPhiE(int i,double u,double v){
  
#ifdef DEBUG
  if(i>nbENP){cerr << "<<<< ERROR >>>> [functionSpace.h] requested phi(i="<< i <<") but nbENP = "<< nbENP <<")" << endl; exit(1);}
#endif
  
  double phi = 0.;
  for(int j=0;j<nbENP;j++) phi += getERefNodesCoef(i,j)*monomeof(u,v,j);
  
  return phi;
  
}
//------------------------------------------------------------------------------
double LagrangeFctSpaceForTriangle::getdPhiE(int i,int dim,double u,double v){

#ifdef DEBUG
  if(i>nbENP)       {cerr << "<<<< ERROR >>>> [functionSpace.h] requested dphi(i="<< i <<") but nbNP = "<< nbENP <<")" << endl; exit(1);}
  if(dim<0 || dim>1){cerr << "<<<< ERROR >>>> [functionSpace.h] requested dphi but dim should be 0(for X) or 1( for Y)" << endl; exit(1);}
#endif
  
  double dphi = 0.;
  for(int j=0;j<nbENP;j++) dphi += getERefNodesCoef(i,j)*dmonomeof(u,v,j,dim);
  
  return dphi;
  
}
//------------------------------------------------------------------------------
double LagrangeFctSpaceForTriangle::monomeof(double u,double v,int i){
  
  return pow(u,pp[i][0])*pow(v,pp[i][1]);
  
}
//------------------------------------------------------------------------------
double LagrangeFctSpaceForTriangle::dmonomeof(double u,double v,int i,int dim){
  
  switch (dim) {
    case 0:
      if(pp[i][0]>0){
        return pp[i][0]*pow(u,pp[i][0]-1)*pow(v,pp[i][1]);
      }else{
        return 0.;
      }
      break;
    case 1:
      if(pp[i][1]>0){
        return pp[i][1]*pow(u,pp[i][0])*pow(v,pp[i][1]-1);
      }else{
        return 0.;
      }     
      break;      
  }
  
}
//------------------------------------------------------------------------------
double LagrangeFctSpaceForTriangle::getPhiF(int i,double u){
  
#ifdef DEBUG
  if(i>nbFNP){cerr << "<<<< ERROR >>>> [functionSpace.h] requested phiF (i="<< i <<") but nbFNP = "<< nbFNP <<")" << endl; exit(1);}
#endif
  
  double phi = 0.;
  for(int j=0;j<nbFNP;j++) phi += getFRefNodesCoef(i,j)*pow(u,(double)j);
  
  return phi;
  
}
//------------------------------------------------------------------------------
void LagrangeFctSpaceForTriangle::setMatrices(){
    
  phiEQP             = DGMatrix(nbEQP,nbENP);
  phiFQP             = DGMatrix(nbFQP,nbFNP);
  Fquv               = DGMatrix(nbFNP,nbFQP);
  Fqduv.push_back(DGMatrix(nbENP,nbEQP));
  Fqduv.push_back(DGMatrix(nbENP,nbEQP));
  dphiUVEQP.push_back(DGMatrix(nbEQP,nbENP));
  dphiUVEQP.push_back(DGMatrix(nbEQP,nbENP));
  
  for(int dim=0;dim<2;dim++){
    vector<DGMatrix> tmp; 
    for(int i=0;i<nbENP;i++) tmp.push_back(DGMatrix(nbENP,nbEQP));
    dphiUVphiEQP.push_back(tmp);
  }
  for(int i=0;i<nbFNP;i++) phiphiFQP.push_back(DGMatrix(nbFNP,nbFQP));

  for(int i=0;i<nbENP;i++){
    for(int q=0;q<nbEQP;q++){
      for(int j=0;j<nbENP;j++){
        for(int dim=0;dim<2;dim++) dphiUVphiEQP[dim][i](j,q) = getdPhiE(j,dim,getEQP(q,0),getEQP(q,1))*getPhiE(i,getEQP(q,0),getEQP(q,1));
      }
    }
  }
  for(int i=0;i<nbFNP;i++){
    for(int q=0;q<nbFQP;q++){
      for(int j=0;j<nbFNP;j++){
        phiphiFQP[i](j,q) = getPhiF(j,getFQP(q))*getPhiF(i,getFQP(q));        
      }
    }
  }

  
  for(int i=0;i<nbEQP;i++){
    for(int j=0;j<nbENP;j++){
      phiEQP(i,j)       = getPhiE(j,getEQP(i,0),getEQP(i,1));
      dphiUVEQP[0](i,j) = getdPhiE(j,0,getEQP(i,0),getEQP(i,1));
      dphiUVEQP[1](i,j) = getdPhiE(j,1,getEQP(i,0),getEQP(i,1));
      Fqduv[0](j,i)     = dphiUVEQP[0](i,j)*EW(i);
      Fqduv[1](j,i)     = dphiUVEQP[1](i,j)*EW(i);
    }
  }
  
  for(int i=0;i<nbFQP;i++){
    for(int j=0;j<nbFNP;j++){
      phiFQP(i,j)       = getPhiF(j,getFQP(i));
      Fquv(j,i)         = phiFQP(i,j)*FW(i);
    }
  }
  
  // --- create Mass matrices --------------------------------------------------
  
  Mass          = DGMatrix(nbENP,nbENP);
  InvMass       = DGMatrix(nbENP,nbENP);
  LumpedMass    = DGMatrix(nbENP,nbENP);
  InvLumpedMass = DGMatrix(nbENP,nbENP);

  for(int i=0;i<nbENP;i++){
    for(int j=0;j<nbENP;j++){
      for(int q=0;q<nbEQP;q++){
        Mass(i,j) += EW(q)*phiEQP(q,i)*phiEQP(q,j);
        if(i==j) LumpedMass(i,j) += EW(q)*phiEQP(q,i);
      }
    }
  }
  
  InvMass.memcpy(Mass);
  InvMass.invert();
  InvLumpedMass.memcpy(LumpedMass);
  InvLumpedMass.invert();
  
}
//------------------------------------------------------------------------------
FunctionSpace* FunctionSpaceFactory::CreateNewFunctionSpace(FctType fctType,int order) {
  
  switch (fctType) {
    case LagrangeTri:{
      return new LagrangeFctSpaceForTriangle(order);
      break;}
    case LagrangeQuad:{
      //NOT YET IMPLEMENTED
      break;}
    default:
      break;
  }
  
}

